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Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations 
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By adding a non-linear core correction to the well established Dual Space Gaussian type pseudopotentials for 
the chemical elements up to the third period, we construct improved pseudopotentials for the Perdew Burke 
Ernzerhof (PBE) [J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)] functional 
and demonstrate that they exhibit excellent accuracy. Our benchmarks for the G2-1 test set show average 
atomization energy errors of only half a kcal/mol. The pseudopotentials also remain highly reliable for high 
pressure phases of crystalline solids. When supplemented by empirical dispersion corrections [S. Grimme, J. 
Comput. Chem. 27, 1787 (2006); S. Grimme, J. Antony, S. Ehrhch, and H. Krieg, J. Chem. Phys. 132, 
154104 (2010)] the average error in the interaction energy between molecules is also about half a kcal/mol. 
The accuracy that can be obtained by these pseudopotentials in combination with a systematic basis set is 
well superior to the accuracy that can be obtained by commonly used medium size Gaussian basis sets in 
all-electron calculations. 



I. INTRODUCTION 

During the last decades, density functional theory 
(DFT) has proven its pivotal role for computational stud- 
ies in the fields of condensed matter physics and quantum 
chemistry. Particularly the Kohn-Sham (KS) formalism 
of DFT has gained enormous popularity as an ab initio 
method applicable to relatively large systems. An essen- 
tial ingredient for many large scale implementations of 
KS-DFT are pseudopotentials which are also frequently 
denoted as effective core potentials. By eliminating the 
strongly bound core electrons pseudopotentials reduce 
the number of occupied electronic orbitals that have to 
be treated in an electronic structure calculation. In addi- 
tion the valence wavefunctions arising from a pseudopo- 
tential are much smoother than the all-electron valence 
wavefunction in the core region, since the orthogonality 
constraints to the rapidly- varying wavefunctions carrying 
core electrons are missing. Since it is not necessary to de- 
scribe rapidly varying wavefunctions the size of the basis 
set used for their representation can be reduced. These 
two factors lead to a significant reduction of the compu- 
tational effort of a pseudopotential calculation compared 
to an all-electron calculation. 

Even though it is well known that the valence electrons 
are responsible for the majority of chemical and physical 
properties of atoms, pseudopotentials have to be con- 
structed very carefully in order to reproduce the proper- 
ties of the all-electron atom accurately. If a pseudoatom, 
i.e an atom described by a pseudopotential, reproduces 
the all-electron behavior accurately for any chemical en- 
vironment the pseudopotential is said to be transferable. 

Pseudopotentials (PSP's) are an essential ingredient 



of most electronic structure codes and different solu- 
tions are implemented in present-day DFT codes. Tra- 
ditional norm-conserving (NC) approaches, e.g.* are for- 
mally the simplest approach, since they give rise to pseu- 
dowavefunctions which lead to a valid charge density. 
By introducing atomic like orbitals as additional basis 
functions any atomic Hamiltonian arising either from an 
all-electron potential or from a norm conserving pseu- 
dopotential^ can be transformed into an Linearized Aug- 
mented Plane Wave (LAPW) like Hamiltonian^"^. The 
widespread Projector-Augmented Wave (PAW) meth- 
ods^ and the ultrasoft pseudopotentials"'^'' are derived by 
such a transformation from an all-electron atom Hamilto- 
nian. The number of required basis functions is reduced 
by this transformation, but the calculation of the charge 
density is more complicated and a generalized eigenvalue 
problem has to be solved even for the case of an orthog- 
onal basis set. For applications in quantum chemistry, 
effective core potentials^"'^"^'^ are often optimized for a 
certain basis set and usually employed for heavier ele- 
ments only. 

In this paper, the Dual Space Gaussian pseudopoten- 
tials of Goedecker-Teter-Hutter (GTH) and Hartwigsen- 
Goedecker-Hutter (HGH)^"* PSP are generalized by the 
inclusion of a Non-Linear Core Correction (NLCC) term. 
These new pseudopotentials are able to provide an ac- 
curacy that is comparable to that of the very best all- 
electron calculations. 

The starting point for understanding why pseudopo- 
tentials work is the subdivision of space into muffin-tin 
spheres centered on the atom in a molecule or solid and 
the remaining interstitial region^^. A non-selfconsistent 
Schrodinger equation can be solved exactly in the inter- 
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stitial region if one knows the scattering properties on the 
surface of the muffin-tin spheres"'^^. The scattering prop- 
erties are typicahy specified by the logarithmic deriva- 
tive as a function of energy. This function is the quo- 
tient of the radial outward derivative and the functional 
value of the wavefunction on the surface of the muffin-tin 
sphere. In this way the boundary conditions are specified 
which are necessary to integrate the Schrodinger equa- 
tion, a second order partial differential equation where 
the amplitude of the solution is fixed by a normaliza- 
tion constraint. A necessary but not sufficient condi- 
tion for a pseudopotcntial to be transferable is therefore 
that the logarithmic derivatives of the all-electron and 
pseudo-atom agree over the relevant energy interval. The 
construction of pseudopotentials is typically done using 
as the reference state a neutral isolated atom which has 
been spherically symmetrized. This symmetrization can 
be achieved by using identical and generally fractional 
occupation numbers for all the orbitals with the same n 
and i quantum numbers, e.g. for the set of 2px, 2py and 
2pz orbitals. The norm conservation conditiou"'^^ ensures 
that the logarithmic derivative function describes well the 
scattering properties of a muffin-tin sphere containing the 
charge distribution of this reference configuration. In a 
selfconsistent calculation the charge distribution changes 
however when the free atom is inserted in a molecule or 
solid and the potential in the muffin-tin region will in gen- 
eral differ from the potential within a muffin-tin sphere 
of the same radius around the reference atom. Hence 
the scattering properties change and the pseudopotcntial 
constructed using the charge distribution in the muffin- 
tin sphere of the isolated atom might not well reproduce 
these modified scattering properties of a new chemical 
environment. Due to screening effects there exists how- 
ever an invariant muffin-tin sphere within which the to- 
tal electronic charge distribution is nearly independent of 
the chemical environment^®. The radius of this invariant 
muffin-tin sphere is a fraction of the covalent bondlength 
and thus considerably smaller than the muffin-tin radii 
used in other methods such as the LAPW method. The 
scattering properties of this invariant muffin-tin sphere 
hardly vary as a function of the chemical environment of 
the atom. If the separable terms of a pseudopotcntial as 
well as the difference between the local part of the pseu- 
dopotentials and the pure coulombic potential remain lo- 
calized within this invariant sphere the pseudopotcntial 
is expected to be highly transferable. This recipe was 
followed in the construction of the GTH^^ and RGE^^ 
pseudopotentials which are indeed well transferable for 
non-spinpolarized systems. 

Despite the fact that the total charge in the invariant 
muffin-tin sphere is nearly identical in different chemical 
environments, the spin polarization is not, as illustrated 
in Fig. 1. Shown are the changes in the radial charge and 
spin densities if one adds half an electron to the unoc- 
cupied spin channel of the 3p-orbital. For spin polarized 
calculations the concept of an invariant muffin-tin sphere 
is therefore not applicable. One possibility to overcome 
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FIG. 1. Difference in the radial spin densities and total charge 
density when adding half an electron to a phosphorus atom 
([A'^e]3s^3p'^'^). The inert core region of the total charge is 
not observed for the individual spin densities. 



this problem is to construct pseudopotentials which have 
an explicit spin dependence^*^ . The other possibility is 
to include nonlinear core corrections (NLCC)^^ in the 
pseudopotcntial. 

In the NLCC schemes the spin and charge densi- 
ties in the muffin-tin sphere are not just the ones from 
the valence electrons treated explicitly by the pseu- 
dopotcntial but they are both the respective sums of 
the valence charge and the core charge given by the 
nonlinear core correction. Since the core electrons 
can be considered to be frozen, i.e to be invariant 
with respect to different chemical environments, this 
core charge is fixed once and for all. It is obvious 
that the spin polarization, i.e. the quantity 0{r) = 
{Pupij) - Pdowniy)) I (Pupi'r) + Pdown{j^)) IS vcry poorly 
represented without a core charge. If for instance all 
valence orbitals are spin up then the spin down charge 
density pdou;n(r) would be zero and the spin polarization 
9 would be equal to one. In the real atom 9 is not equal 
to one in the core region since the core electrons are never 
spin polarized. Since the density of the core electrons is 
much larger than the density of the valence electrons in 
the core region the spin polarization actually has typi- 
cally small values. These correct small values of 9 are 
reestablished by the core charge of a NLCC pseudopo- 
tcntial and exchange correlation functionals can provide 
reliable total energies. Nonlinear core corrections have 
therefore the potential to substantially improve the de- 
scription of spin polarized states. 

Whereas previous implementations of NLCC pseu- 
dopotentials tried to faithfully represent the core 
charge, we follow here another approach. In the spirit of 
the GTH pseudopotentials where all terms have simple 
parametrized analytical forms we also represent the core 
charge density just as one single Gaussian. The ampli- 
tude and width of this Gaussian core charge distribution 
are then optimized by a fitting procedure in the same 
way as the other parameters of the pseudopotential. 
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II. METHODOLOGY 

The procedure for constructing the NLCC pseudopo- 
tcntials is very similar to the one used for the construc- 
tion of the GTH and HGH pseudopotential. In con- 
trast to the original GTH and HGH pseudopotentials 
which were fitted to a single atomic configuration, the 
new NLCC pseudopotentials are fitted not only to the 
ground state but also to several excited and ionized elec- 
tronic configurations where half an electron is added or 
removed possibly to or from different valence orbitals. 
The atoms are always considered to be spherically sym- 
metric, but some of the configurations used for the fit 
have a spin polarization. The parameters of the dual 
space Gaussian ansatz^^ are fitted such that: 

• The occupied and first few unoccupied valence 
eigenvalues of the all-electron and pseudo atom 
match for all configurations used in the fit. 

• The charge inside the inert region of the pseudo 
atom matches the all-electron valence charge in the 
same region for all the orbitals for which the eigen- 
values are matched for all configurations used in the 
fit. This means that the pseudopotential is norm 
conserving for all the configurations used in the fit. 

• A high precision of 10~^ a.u. is achieved for va- 
lence eigenvalues and charge integrals of the non- 
polarized ground state, whereas only a moderate 
precision of lO"** a.u. is enforced for all other or- 
bitals and configurations considered. 

• The total energy differences of the all-electron atom 
are reproduced for all configurations used in the fit. 

• The spin polarization energies of the all-electron 
atom are reproduced for all spin polarized configu- 
rations used in the fit. 

Since the considered quantities are fitted for several con- 
figurations atomic transferability is already built in to 
these new pseudopotentials by construction. 

The core charge is represented by a single Gaussian 
which is optimal for numerical efficiency. It is initialized 
such that it approximates well the physical core charge 
density and it is held fixed during an initial stage of the 
fit. Then both the amplitude and the width of the Gaus- 
sian are released and considered as fitting parameters. 
As a consequence the total amount of core charge and 
the width can differ from the physical value. 

The parameters of the core charge constitute thus a 
small set of only two additional degrees of freedom. Yet 
this allows to optimally reproduce atomic polarization 
energies without degrading the transferability and accu- 
racy of other atomic properties. It was found that the 
inclusion of a more complicated core charge is not benefi- 
cial. Furthermore, it should be emphasized that the novel 
NLCC pseudopotentials are not harder than their HGH 
counterparts. The smoothness of the core charge seems 



to play an important role for the fact that the hardness 
is not affected, and roughly the same grid spacings or 
energy cutoffs can be employed as for conventional HGH 
pseudopotentials. 

In particular, pseudopotentials with NLCC were gen- 
erated for boron, carbon, nitrogen, oxygen, fiuorine, alu- 
minium, silicon, phosphorous, sulfur and chlorine. Very 
weak spin dependences are expected for the rare gasses, 
and for all remaining chemical elements up to the third 
row, NLCC are found to be unnecessary, as HGH pseu- 
dopotentials are available that either include semicores 
(sodium and magnesium) or leave no core states at all 
(hydrogen, lithium and beryllium). For the special case 
of hydrogen, it was found that the multi configuration 
fit gave slightly improved results even though obviously 
no core charge was added. Since the focus of this paper 
is on systems made out of light elements, no relativistic 
effects such as spin-orbit coupling were included in the 
pseudopotentials. 



III. RESULTS 

To assess the accTiracy of the new pseudopotentials ex- 
tensive calculations were performed for different test sets. 
The accuracy of covalent bond formation energies was 
examined for the standard G2-1 test set^'^"^^. For the as- 
sessment of the accuracy of non-bonded interactions the 
8222'^ test set was used. To check the performance for 
materials under high pressure we chose carbon, silicon, 
silicon carbide and boron nitride as test systems. 

All pseudopotential calculations were done with the 
BigDFT package^^. The BigDFT code uses a system- 
atic wavelet basis set which allows to obtain the exact 
density functional solution with arbitrarily small error 
bounds. The parameter were set such that an accuracy 
of at least 10~® Hartree was obtained. The LibXC li- 
brary^^ is used within BigDFT for the evaluation of the 
exchange correlation functional. Semi-empirical van-der 
Waals corrections were added in BigDFT according to 
the DFT-D22 and DFT-DS^ methods for the calculations 
of the S22 test set. 

To obtain reliable all-electron reference values for 
the atomization energies of the G2-1 test set, we per- 
formed all-electron calculations with the NWChem soft- 
ware package^*' using one of the largest available Gaus- 
sian type basis sets, namely an augmented correlation 
consistent polarized valence quintuple zeta Gaussian type 
basis set (aug-cc-pV5Z). Care was taken to disable sym- 
metry detection and to check for the lowest energy spin 
multiplicity. For the chemical elements Li, Be Na and 
Mg, the aug-cc-pV5Z set was not available, so the cor- 
responding quadruple zeta set (aug-cc-pVQZ) was used 
to compute the atomization energies of Li2, LiF, BeH, 
Na2, NaCl, MgH and Mg2. To obtain the atomization 
energies of the relaxed molecules, geometry optimizations 
were carried out using the very same basis set. 

Atomic all-electron calculations of the spin polariza- 
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tion energies were done with our non-spherical atomic 
code, which expresses the wavefunctions as a product 
of spherical harmonics and radial functions. The radial 
function are given numerically on a logarithmic grid. The 
settings were chosen such that a precision of at least 10~® 
Hartree can be obtained for the total energy. This re- 
quired angular integration grids of 232 points and multi- 
pole representations up to ^ = 4. The atomic LSDA ref- 
erence energies from the National centre of science and 
technology (NIST)"^^ where reproduced within the given 
precision for all elements considered. 

Wc calculated the atomization energies also with the 
three different sets of PAW^ potentials available in 
VASP^^. Those PAW potentials are derived from the 
all-electron atomic Hamiltonian and aim at all-electron 
accuracy. In order to obtain the required high precision 
some parameters had to be set to tighter values than the 
default values. We had to use for the general accuracy 
(PREC = High Accurate and lasph = true) to ac- 
tivate nonsphcrical gradient corrections inside the PAW 
spheres. It was carefully checked that the calculations 
were converged with respect to the size of the periodic 
simulation cell. Furthermore, care was taken that the 
correct spin multiplicity and non-fractional occupations 
were produced. Hard PAW potentials were available for 
all required elements except for Li, Be, Na and Mg, for 
which scmicore potentials were used instead. For com- 
parison, all energies wcirc; recomputed with a set of default 
potentials. The third set consists of soft potentials for the 
elements B, C, N, O and F and default potentials oth- 
erwise. For the periodic solids, all-electron calculations 
have been performed using the full-potential linearized 
augmented plane wave (FLAPW) and augmented plane 
wave plus local orbitals (APW + lo) methods as imple- 
mented in the WIEN2k'^'^ software package. We used a 
reduced mufRn-tin radii for all atomic sorts in order to 
avoid their overlap up to the highest studied pressures. 
The sphere radii were kept fixed throughout the whole 
set of lattice parameters to obtain the best possible er- 
ror cancellation. Semicore states were treated as valence, 
because high compressions can lead to an overlap of their 
wavefunctions, which will give a contribution to the en- 
ergy. Inside the spheres, the partial waves were expanded 
up to LMAX = 10. The number of plane waves was lim- 
ited by a cutoff parameter RMTKmax = 9.0 for all the 
compounds under consideration. The charge density was 
Fourier expanded with Gmax = 14 a.u. For the ma- 
jority of the systems wc used a very dense k-points grid 
(15x15x15) to ensure total energy convergence. 

All the calculations were done at zero electronic tem- 
perature, i.e. no Fermi smearing was used. Zero point 
energies were not included in any of our results. 



A. Atomization energies of the G2-1 test 

Atomization energies are frequently used to assess the 
quality of various exchange correlation functionals as well 



as other approximations used in electronic structure cal- 
culations. The Gaussian G2-1 test set^'^"'^® is a standard 
benchmark set of 55 molecules in this context. Since 
this test set does not contain molecules with the chemi- 
cal elements B, Al and Mg we added the molecules BH, 
BH2, AlH, A1H2, Mg2 und MgH. We used this aug- 
mented test set to compare our pseudopotential results 
with all-electron calculations. Because of Hund's rule 
most isolated atoms are strongly spin polarized. When 
an atom is inserted into a molecule or solid, its spin po- 
larization is typically strongly reduced. Since standard 
pseudopotentials arc based on a non-spin polarized ref- 
erence configuration they can typically better describe 
atoms in molecules or solids than isolated atoms them- 
selves. Since the atomization energy is the difference be- 
tween the total energy of the molecule and the sum of 
the total energies of its constituent isolated atoms, the 
largest contribution to the error in the atomization en- 
ergy of a pseudopotential calculation comes actually from 
the atomic energies. 

The atomization energies of the molecules in the G2- 
1 test set were first computed using conventional HGH 
pseudopotentials^^ for the PBE exchange correlation 
functional. A comparison with all-electron data is shown 
in figure 2. The spin multiplicity of systems with a net 
magnetic moment arc indicated in brackets and omitted 
for closed shell systems. Deviations of ±1 kcal/mol are 
indicated with a (green) shading to relate the errors to 
the requirements for chemical accuracy. 

It is found that the direct computation of the electronic 
atomization energies with the conventional pseudopoten- 
tials leads to significant disagreement with the results 
obtained in an all-electron calculation. An rather high 
mean absolute deviation (MAD) of 6.83 kcal/mol to the 
electron reference values for all 55 molecules in the G2-1 
set is found. However, the main contribution to the error 
in the atomization energies comes from the estimation of 
the energy of the isolated atoms. Therefore, the atomiza- 
tion energies can be improved significantly by a two step 
procedure where the atomization energies are calculated 
as a sum of two terms. The first term is the energy differ- 
ence between the molecule and the sum of the total en- 
ergies of isolated, spherical and non-spinpolarizcd atoms. 
It thus can be considered as the atomization energy with 
respect to a set of non-physical atoms. This energy dif- 
ference is calculated with the HGH pseudopotentials and 
is fairly accurate since no strong spin polarizations are 
involved. The second term is the difference in total en- 
ergy between the real, i.e non-spherical and spin polar- 
ized, atom and the previously defined non-physical atom. 
This second term is calculated with our all-electron pro- 
gram for non-spherically symmetric atoms and is there- 
fore exact. Since the atomic spin polarization energies 
and energy terms for breaking the spherical symmetry 
are only a property of the atoms they can be considered 
as a set of atomic correction terms for the accurate cal- 
culation of atomization energies. The atomic correction 
terms for the chemical elements considered in this study 
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TABLE L Atomic correction terms in kcal/mol as used for 
the two step procedure. 



are listed in Table I. 

It has to be stressed that these atomic spin polarization 
energies drop out in most instances such as in the calcu- 
lation of energy differences in a chemical reaction where 
only molecules are involved. Using this two step scheme, 
the errors in the atomization energies are decreased con- 
siderably to a MAD of 1.56 kcal/mol. Because of the can- 
cellation effect, this is the accuracy that can be expected 
in the majority of energy differences calculated with the 
standard HGH type pseudopotentials. The above MAD 
value was obtained with the bond lengths and angles 
fixed as given on the computational chemistry compari- 
son and benchmark database'^^. Using instead the equi- 
librium geometries of each method, i.e. the HGH pseu- 
dopotential and the all-electron calculation, the MAD 
value is slightly decreased to 1.52 kcal/mol. This last 
value is actually more relevant in practice since atomiza- 
tion energies for an unknown system necessarily have be 
calculated with theoretically determined geometries. 

The new Gaussian type pseudopotential with a NLCC 
can however still considerably improve the accuracy with- 
out the need of using a two step procedure. The error of a 
direct computation of the atomization energies decreases 
to a MAD of only 0.52 kcal/mol. Using the equilibrium 
geometries obtained with the new NLCC pseudopoten- 
tial the error drops again slightly to 0.51 kcal/mol. More 
important than this small improvement of the MAD is 
the fact that the result could be improved for the few 
molecules where the error was well above the average. 

In figure 3 the accuracy of the HGH pseudopotentials 
with NLCC is shown for relaxed molecular geometries 
and compared with the results of PAW calculations pub- 
lished by Paier et. al."^^. In this work hard PAW po- 
tentials were used and we were able to reproduce their 
results. In essence, the absolute errors of the new NLCC 
pseudopotentials are comparable with those using hard 
PAW potentials. As shown in Figure 4 the accuracy how- 
ever goes down significantly when one uses the default or 
even the soft PAW potentials of the VASP package^^. 
Furthermore, the same figures shows the discrepancies 
between all-electron results obtained in two large Gaus- 
sian basis sets while keeping the molecular geometries 
fixed. Even at this size the differences between the two 
basis sets are not negligible compared to the deviations 
to other methods and the accuracy of the pseudopoten- 
tial method is indeed close to the discrepancies between 
different choices of all-electron reference values. This is 
quite surprising given the fact that these simple chemi- 
cal compounds show only straightforward covalent type 
bonding properties which are certainly easier to describe 
with a Gaussian basis set than other more complex bond- 
ing patterns. It has also to be stressed that the com- 



FIG. 2. Accuracy of PBE^ atomization energies computed 5 
with HGH pseudopotentials. Explanations are in the text. 
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TABLE II. Deviation measures of the electronic atomization 
energies in kcal/mol for the 55 molecules of the G2-1 test set 
compared to the all-electron result obtained in the aug-cc- 
pV5Z basis set. All geometries are fixed. 



putational cost rises very steeply when one goes from a 
medium size basis set to these very large basis sets. This 
is in contrast to the wavelet method where a modest de- 
crease of about 15 percent in the grid spacing h results 
in an gain of a factor of ten in accuracy because of the 
high order convergence rate of /i^^. 

The accuracy problems of Gaussian basis sets become 
even more evident if one employs medium size or small 
standard basis sets in an all-electron calculation. The 
6-31G, 6-31++G*, 6-31+G** and 6-311-F-HG(3df,3pd) 
basis sets were employed to compare the relative accuracy 
of the pseudopotential method with the incompleteness 
of and disagreement between standard Gaussian basis of 
various sizes. Figure 4 clearly shows that the accuracy 
obtained with these basis set is considerably lower than 
the accuracy with the NLCC pseudopotentials or also 
with the standard HGH pseudopotential within the two 
step procedure described above. 

A summary of the deviations in the atomization ener- 
gies averaged over the molecules of the G2-1 test set is 
given for fixed and relaxed geometries in tables II and 
HI, respectively. Indicated are the MAD, RMSD, mean 
signed deviation (MSD), maximum absolute deviation 
(maxAD) and minimum absolute deviation (minAD). 

The last row of table HI describes the change in the 
all-electron reference values when going from the fixed, 
experimental (CCCBDB) to relaxed geometries in the 
aug-cc-pV5Z basis set. This gain in energy upon geom- 
etry relaxation is significant compared to the assessed 
accuracy of the pseudopotential based methods, which 
are found to be very reliable for geometry optimizations. 



B. Accuracy of the equilibrium geometries 

In order to compare the accuracy of the equilibrium ge- 
ometries of the pseudopotential and all-electron calcula- 
tions, the optimized geometry of each molecule is aligned 
with its all-electron counterpart, such that the RMSD is 
minimized^^. The resulting RMSD values are shown in 
figure 5. It is observed that conventional HGH pseu- 
dopotentials yield already very good agreement with the 
all-electron data. Nevertheless, the inclusion of NLCC 



FIG. 3. Comparison of PBE atomization energies from 6 
NLCC-HGH pseudopotentials with other methods. Expla- 
nations are in the text. 




FIG. 4. Comparison of PBE atomization energies from 7 FIG. 5. RMSD values of the experimental and relaxed ge- 
NLCC- HGH pseudopotentials with less accurate methods. ometries with respect to those relaxed with all-electron cal- 
Explanations are in the text. culations in the aug-cc-PV5Z set. 
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TABLE III. Deviation measures of the electronic atomization 
energies in kcal/mol for the 55 molecules of the G2-1 test set, 
where all molecular geometries are optimized for each method 
considered. For comparison, PAW data are extracted from 
work of Paier et. al^®. The last row gives the change of the 
all-electron energy upon geometry relaxation. 
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leads to a systematic improvement of the equilibrium ge- 
ometries. 

It can be noted that in a previously mentioned work'^^, 
a similar test was carried out for the bond lengths of 
some dimers in order to verify the accuracy of the PAW 
method. It is found that our pseudopotential approach 
yields geometry data of at least the same or even better 
accuracy, and that the high precision is maintained when 
moving to more complicated geometries. 



C. Evaluation of pressure of extended systems 
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Next we benchmark pseudopotential (PSP) calcula- 
tions for extended systems. A few crystalline systems 
made of light elements (diamond carbon, Silicon Car- 
bide, Bulk Silicon and Boron Nitride) were selected and 
the pressure at a given lattice parameter was then com- 
pared between different approaches. 

The details on how the stress energy tensor can be 
calculated in GGA for NLCC terms are given in the ap- 
pendix. Figure 6 shows the difference between the LAPW 
results and the PSP results at the same lattice param- 
eter. In addition we also show the results for the hard 
PAW potentials. In order to show the relative accuracy 
of the pressure we specify the lattice constant along the 
X axis of the figure in terms of the pressure. At the high- 
est pressures the lattice constant is reduced by about 10 
percent compared to the lattice constant with zero pres- 
sure. With NLCC PSP, an excellent relative accuracy 
of about 10^"^ is found. In this case, it can be noticed 
that the inclusion of a NLCC term improves further the 
results even though the systems under pressure are not 
spin-polarized. Results of similar quality can be obtained 
within the hard PAW scheme described above. 



D. Dispersion-corrected functionals 

Long range van der Waals interactions are missing in 
all standard LDA and GGA density functionals. Adding 
semiempirical classical van der Waals interactions has 
however recently been demonstrated to give a rather ac- 
curate description of weakly bonded systems and is now 
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FIG. 6. Comparison of pressures. 
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MAD 


RMSD 


MSD 


maxAD 


minAD 


NLCC-no corr. 


2.59 


3.61 


2.59 


10.08 


0.05 


HGH Krack-no corr. 


2.64 


3.66 


2.64 


10.18 


0.01 


NLCC-D2 


0.51 


0.64 


-0.39 


1.57 


0.05 


HGH Krack-D2 


0.50 


0.58 


-0.34 


1.25 


0.13 


NLCC-D3 


0.48 


0.64 


-0.03 


1.14 


0.05 


HGH Krack-D3 


0.47 


0.64 


0.02 


1.95 


0.03 


aug-cc-pVDZ-D2 


1.05 


1.15 


-1.05 


2.33 


0.49 


aug-cc-pVTZ-D2 


0.53 


0.68 


-0.51 


1.60 


0.02 


aug-cc-pVTZ-D3 


0.44 


0.57 


-0.14 


1.43 


0.07 



TABLE IV. Deviation measures in kcal/mol for the S22 test 
set with respect to CCSD(T) calculations. For the PBE XC 
functional, PSP and all-electron calculations are compared 
including semiempirical dispersion corrections (D2 and D3). 



frequently used. We will therefore examine the accuracy 
of the semiempirical models in the context of our pseu- 
dopotential calculations with a systematic wavelet basis 
set. 

In BigDFT, we implemented two semiempirical models 
to correct dispersion energies and energy gradients DFT- 
D2^ and DFT-D3'^. The parameters of these models were 
separately fitted for each exchange correlation functional 
based on thermochemical data for weakly interacting sys- 
tems. Since BigDFT uses a wavelet basis and pseudopo- 
tentials Figure 7 and Table IV show the comparison of in- 
teraction energies of the benchmark database 822^'', with 
a reference calculation using Coupled Cluster CCSD(T) 
in the complete basis limit (CBS)"^*. The inclusion of 
dispersion correction D2 into BigDFT clearly improves 
the description of weak interactions within PBE, even 
though the S22 data set was not used as the fitting data 
set. The root-mean-square-deviation (RMSD) between 
the CCSD(T) reference values and the NLCC-DFT inter- 
action energies is 0.58 kcal/mol. The absolute maximum 
difference corresponds to acetic acid dimer (C00H)2, 
where the overestimation is 1.57 kcal/mol, that means 
an 8% of the total interaction energy. The largest rela- 
tive error of 35 % is found for the methane dimer whose 
interaction energy is only 0.2 kcal/mol. The errors for 
these systems are comparable to those that are obtained 
when PBE-D2 and PBE-D3 are used with a large ba- 
sis set (aug-cc-pVDZ and aug-cc-pVTZ). On average the 
PBE-D2 scheme performs better with BigDFT than with 
any Gaussian basis set, while the PBE-D3/BigDFT re- 
sults are comparable to the results obtained with PBE- 
D3/aug-cc-pVTZ electron calculations. 



IV. DISCUSSION AND CONCLUSIONS 

We have shown that our new NLCC PSP's give very 
high accuracy for a wide range of applications. In par- 
ticular they give atomization energies with chemical ac- 
curacy compared to all-electron calculations for the G2-1 
test set. This accuracy can easily be obtained with a 
systematic basis set such as wavelets where one has to 




I I I I I I I I 

2 1.5 1 0.5 0-0.5 -1 -1.5 -2-2.5 
deviation to CCSD(T)/CBS in kCal/mol 



FIG. 7. Comparison of S22 test set results between PSP and 
all-electron calculations within PBE XC functional. 



change only a single parameter to obtain arbitrarily high 
accuracy. Obtaining such a high accuracy with Gaussian 
basis sets requires using the largest available basis sets 
and is therefore frequently not feasible in practice. Con- 
trary to a widespread belief, PAW calculations do not 
necessarily give all-electron accuracy. Soft PAW poten- 
tials can actually lead to appreciable errors. Well con- 
structed hard PAW potentials on the other hand give 
very high accuracy and are together with our new norm- 
conserving pseudopotentials in practice the only feasible 
way to highest quality results for large systems. 
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V. APPENDIX: NLCC HGH 
PSEUDOPOTENTIALS IN KOHN-SHAM DFT 
FORMALISM 



The PSP format is based on HGH-Krack form^^: 

VpSP = Vloc + Vnl , 

where the first part is a local potential 



(1) 



Viocir) 



+ exp 




2k-2\ 



(2) 



and the second part the non-local term which is sepa- 
rated into different channels Vni = Ve{r, r'), each one 
defined in terms of separable projectors 



n<2 e 

Ve{r,r')=J2 E pf^W /^f,- pf (r') 

i,j=l m=—i 



V2 r^^+'e 



The core charge pc of the new PSP's is given by 

7—7 



Pc{r) = Cc, 



(\/27rrcore)' 



(3) 
(4) 

(5) 



The pseudopotentials parameters according to equa- 
tions (2) to (5) are given in table V. 

The core density is then used in the Kohn-Sham total 
energy expression as follows: 

i 

- Eh \p] + i^xc [p + j drp(r)14e [p + Pc] (r) , (6) 

where E^c and 14c[?i] = ^^sn"^^ energy and 

potential respectively, Vh is the Hartree potential and 
the tpi^s are KS wavcfunctions, whose summed squares 
give the valence density p = IV'iP- Eq-6 ensures 
Hellmann-Feynman condition at self-consistency, ^S<s = 
0. 

The contribution to the stress tensor coming from 
the XC term with NLCC can be shown to be given by 



OT^I = 5o,i3E^^[p + Pc] - (5a/3 J drV^cb + Pc](r)p(r) 

-I- J drV:^c[p + Pc]{r)r(ad0)Pc{r) 
-/dr(n(rMn](^)(r)^)a,)n(r) 



(7) 
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18.37181 
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0.26844 


0.10004 




rp 








0.25234 


0.44314 






C-core 




F 


9 


7 




Z 








0.20610 
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2.79309 
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C2 




0.19518 


23.47047 
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0.00000 
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0.48775 


0.38780 
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14 
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0.33000 


-0.07846 
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rioc 


Cl 
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2.87392 
0.00000 


0.02559 
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rs 
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0.48800 


2.47963 
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0.44279 
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15 
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0.34000 


-1.62258 


-0.72412 


rioc 


Cl 


C2 




0.38209 


3.47754 


-0.01267 
3.47461 


rs 
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hh 




0.43411 


3.37859 




rp 
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TABLE V. Pscudopotcntial parameters of HGH potentials 
with NLCC. The ionic charge Zion, local radius rioc and coef- 
ficients Ck define the local part (2), while the separable part 
(3) is determined by the localization radii r^ and coefficients 
hij. The parameters for the core charge (5) 3ir6 Ccore and 

rcore . 



n=p+Pc 



10 



where is the supercell vohimc and £[n]*-^^ = 
de[n]/d{\S/n\). The formula shows that the gradient of Pc 
is needed to evaluate T^p, even for a LDA computation. 
A detailed derivation of DFT of stress (without NLCC) 
was shown by Dal Corso and Cresta^^. 
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